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Abstract 


An  efficient  method  for  unstructured  mesh-to-mesh  interpolation  is  described.  This 
method  uses  a  binary  space  partitioning  tree  to  sort  the  elements  of  the  source  mesh.  Using  this 
tree  data  structure  the  source  mesh  elements  ean  be  efficiently  searched  to  find  the  nearest 
clement  for  each  of  the  destination  mesh  points.  Once  found,  the  nearest  source  mesh  element  is 
used  to  compute  baryeentric  coordinates  which  are  then  used  as  weighting  coefficients  for  the 
interpolation. 


Administrative  Information 

The  work  described  in  this  report  was  performed  by  the  Computational  Hydromechanics 
Division  (Code  5700)  of  the  Hydromechanics  Department  at  the  Naval  Surface  Warfare  Center, 
Cardcrock  Division  (NSWCCD).  The  work  was  funded  by  the  Office  of  Naval  Research 
(Program  Officer:  Dr.  Ki-Han  Kim)  under  the  Force  Protection  Applied  Research  Program  (PE 
602123N). 


Background 

This  report  describes  the  code  and  algorithms  used  to  perform  a  one  way  coupling  of  a 
Computational  Fluid  Dynamics  (CFD)  solver,  and  a  Computational  Structural  Dynamics  (CSD) 
solver.  This  code  was  developed  at  the  Naval  Surface  Warfare  Center,  Carderock  Division  under 
the  FY08  ONR  funded  Crashbaek  program.  Specifically  the  methods  used  to  transpose  the  time 
varying  fluid  dynamic  surface  pressure  field  computed  by  the  CFD  solver,  to  the  wetted  faces  of 
the  finite  element  mesh  used  in  the  CSD  solver  are  described.  In  this  report  the  CFD  surface 
mesh  is  referred  to  as  the  source  mesh,  and  the  CSD  mesh  as  the  destination  mesh.  Transposing 
the  surface  pressure  field  requires  that  for  each  point  in  the  destination  mesh,  a  nearest  element  is 
found  in  the  source  mesh  from  which  all  desired  scalar,  and  vector  quantities  may  be 
interpolated. 

Binary  Space  Partitioning  Tree 

Due  to  the  relatively  large  mesh  sizes  involved  (typically  greater  than  1x10^)  a  Binary 
Space  Partition  (BSP)  tree,  a  type  of  binary  search  tree,  is  used  to  sort  the  elements  of  the  source 
mesh.  This  offers  an  efficient  mechanism  for  identifying  the  nearest  clement  in  the  source  mesh 
to  a  given  point  in  the  destination  mesh.  Given  the  large  mesh  sizes  and  that  each  destination 
mesh  point  needs  an  interpolated  value,  the  time  invested  in  building  the  BSP  tree  is  considered 
acceptable.  The  effort  involved  in  building  the  BSP  tree  is  expected  to  be  on  the  order  of  nlog(n) 
where  n  is  the  number  of  elements  in  the  source  mesh.  Once  the  tree  is  built  the  time  necessary  to 
search  for  the  nearest  clement  is  expected  to  be  on  the  order  of  log(n)  [1]. 

Building  the  Tree 

For  this  application  a  variant  of  a  BSP  tree  known  as  a  generalized  pseudo  kd-tree  [2]  is 
used.  To  build  the  tree  the  source  mesh  elements  arc  recursively  divided  into  two  bounding 
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boxes.  A  bounding  box  is  a  rectilinear  box  aligned  with  the  coordinate  axis.  It  represents  the 
minimum  bounding  volume  of  all  elements  associated  with  the  node  of  the  tree.  It  could  be 
defined  in  terms  of  Xmin,  Xmax,  Ymin,  Ymax,  and  Zmin,  Zmax,  howcvcr  for  simplicity  the  vectors 
Pmin(Xn,in,Yn,in,  Z^jn),  and  Pniax(Xmax,Ymax,Zn,ax)  arc  defined.  Each  division  is  based  on  cutting  the 
parent  node’s  bounding  box  across  its  longest  axis,  and  through  the  average  value  of  all  the 
element  centers.  Using  the  element  centers  simplifies  the  division,  and  sorting  of  the  elements. 
Each  node  of  the  tree  contains; 

•  A  pointer  to  an  array  of  all  elements  {^Elements) 

•  An  index  to  the  first  element  in  the  tree’s  node  (/„) 

•  The  number  of  elements  in  the  node  (NumberOJSubElements) 

•  Two  vectors  which  represent  a  bounding  box  (Pm/n,  and  Pmar) 

•  Pointers  to  both  sub  nodes  of  the  tree  (5'(j,  S/) 

The  original  array  of  all  elements  is  used  and  never  copied  during  the  process  of  building  the 
tree.  Each  sub-node  simply  points  to  a  sub-section  of  the  array.  Elements  within  the  array  arc 
than  reorganized/sorted  by  the  sub-nodes.  The  root  node  of  the  tree  is  created  by: 

•  Loading  all  elements  from  the  source  mesh  into  the  element  array 

•  Setting  lo  to  zero 

•  Setting  NumberOfSubElements  to  the  total  number  of  elements  in  the  source  mesh 

•  Looping  through  all  nodes,  of  each  element,  and  calculating  Pmin,  and  P„ax 

•  Setting  the  Pointers  So,  and  Si  to  Null 

The  tree  is  grown  by  recursively  subdividing  the  root  until  the  NumberOfSubElements  is  50  or 
less.  Nodes  which  have  50  or  less  elements  become  terminal,  or  leaf  nodes.  The  procedure  for 
recursively  subdividing  the  root  node  is: 

If  and  only  if  the  NumberOJSubElements  is  greater  than  50: 

•  Find  the  longest  edge  of  the  bounding  box 

•  Find  the  average  of  all  the  element  centers 

•  Divide  the  node’s  bounding  box  with  a  plane  that  cuts  orthogonally  through  the 
longest  edge  of  the  bounding  box  and  through  the  average  of  the  element  centers 

•  Sort  the  elements  in  the  current  cell  so  that  all  elements  whose  centers  arc  below  the 
plane  go  in  front  of  all  the  elements  that  arc  not 

•  Create  new  So' 

O  So  ^lo  In 

o  So—*NumberOJSubElements  =  Number  of  elements  found  below  the  plane 
o  Loop  through  all  nodes,  of  each  element  in  So,  to  calculate  So—*Pmin,  and 

So  *P max 

•  Create  new  Si'. 
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o  S/  — =  So—*NumberOJElements 

o  Si—>’NumberOfSubElements='Numhcr  of  elements  not  found  to  be  below  the 
plane 

o  Loop  through  all  nodes,  of  eaeh  element  in  Si,  to  caleulate  the  Si—*P„,i„,  and 

Si  *P max 

•  Repeat  proeedure  for  So,  and  Si 

Searching  for  the  Nearest  Element 

Onee  the  tree  is  populated  it  is  used  to  seareh  for  the  nearest  element  in  the  souree  mesh  to 
a  destination  point.  Given  the  eoordinates  of  a  destination  point,  the  distanee  between  it  and 
each  sub-node’s  bounding  box  is  calculated,  (zero  being  returned  if  the  point  is  inside  the 
bounding  box).  If  one  has  a  eloser  distanee  it  is  searehed  first,  otherwise  Si  is  arbitrarily 
searehed  first.  This  proeess  is  repeated  reeursively  until  a  terminating  node  is  reaehed;  at  this 
point  eaeh  element  assoeiated  with  this  terminating  node  is  tested  individually  to  find  whieh  is 
elosest  to  the  destination  point.  The  index  of  the  elosest  element,  and  the  distanee  that  it  is  from 
the  destination  point  is  retained  during  the  seareh  proeess  to  eliminate  braneh  searehes  whose 
bounding  box  is  farther  from  the  destination  point  than  the  eurrent  nearest  element.  For  our 
partieular  applieation  the  souree  mesh  eontained  only  triangular  elements. 

Finding  the  Distance  to  a  Triangle 

Before  aeeurately  eomputing  the  distanee  to  a  souree  triangle  two  preliminary  tests  are 
performed  in  an  attempt  to  minimize  the  number  of  triangles  that  are  rigorously  tested,  thus 
saving  additional  expensive  ealeulations.  First  eaeh  destination  point  will  have  an  average 
normal  ealeulated  for  it  based  on  all  the  destination  mesh  elements  that  it  is  part  of  Each 
triangle  is  tested  to  find  out  if  the  dot  produet  of  the  average  nonnal  of  the  destination  point  and 
the  normal  of  the  souree  triangle  is  greater  than  zero.  Considering  only  souree  elements  for 
whieh  this  dot  produet  is  greater  than  zero  ensures  that  only  the  souree  elements  that  are  faeed  in 
the  same  general  direetion  as  the  surfaee  of  the  destination  point  arc  used  for  interpolation.  The 
second  test  is  to  sec  if  its  distance  to  the  plane  of  the  triangle  is  greater  than  the  distanee  to  the 
closest  triangle  so  far.  If  this  distance  is  less  than  the  distance  to  the  current  nearest  triangle  then 
more  rigorous  ealeulations  arc  required  to  determine  the  actual  distanee  to  the  triangle. 

To  begin  calculating  the  actual  distance  to  the  triangle  a  test  is  performed  to  find  out  if  the 
destination  point  is  over  the  triangle.  This  test  is  performed  by  considering  a  tetrahedron  with 
the  destination  point  (Pj)  at  the  peak,  and  the  points  to  the  source  triangle  (Po,  P|,  P2)  at  the  base. 
If  the  point  Pd  projects  into  the  base  then  the  distance  to  the  plane  of  the  triangle  previously 
calculated  is  in  fact  the  distanee  to  the  triangle.  If  Pd  projects  outside  the  base  then  one  of  the 
faces  will  be  obtuse  to  the  base,  and  the  distance  from  the  edge  to  Pd  is  required.  Specifically  if 
the  normal  defined  by  souree  triangle  (Po,  Pi,  P2)  dotted  with  the  normal  defined  by  triangle  (Pd, 
Po,  Pi)  is  less  than  zero  then  that  face  of  the  tetrahedron  is  obtuse  to  the  base  and  another 
calculation  will  be  required  to  find  the  distance  to  the  Po,  P|  edge.  Similarly,  if  the  nonnal 
defined  by  triangle  (Po,  Pi,  P2)  dotted  with  the  normal  defined  by  triangle  (Pd,  P|,  P2)  is  less  than 
zero  then  another  calculation  will  be  required  to  find  the  distance  to  the  P|,  P2  edge.  Also,  if  the 
normal  defined  by  triangle  (Po,  Pi,  P2)  dotted  with  the  normal  defined  by  triangle  (Pd,  P2,  Po)  is 
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less  then  zero  then  the  distanee  to  the  P2,  Po  edge  is  required.  Figure  1  eontains  pseudo  eode  for 
this  method.  This  pseudo  eode  assumes  that  VeetorType  defines  the  dot  produet  using  the 
operator  and  the  eross  produet  using  the  “1”  operator. 

Finding  the  Distance  to  an  Edge 

If  during  the  proeess  of  ealeulating  the  distanee  to  a  souree  mesh  triangle  it  is  found  that 
the  destination  point  eannot  be  projeeted  onto  the  triangle  then  the  distanee  to  the  elosest  edge  is 
eomputed.  For  example,  assuming  that  the  edge  Po,  P|  is  the  elosest  edge,  the  first  step  is  to  test 


Float 64  BinaryTreeType : : 

DistanceFromTheTriangle (const  VectorType& 

const  VectorType& 
const  ElementTypeS 
SearchResultsType& 

{ 

VectorType&  PO  =  Element. PO 
VectorType&  Pi  =  Element, PI 
VectorType&  P2  =  Element, P2 


Pd, 

PdNormal , 

Element, 

Results) 


VeetorType  BaseNormal  =  NormalToTriangle ( PO,  PI ,  P2 ) ; 

Float64  d=Abs (BaseNormal* (Pd-PO) ) ;  //  d=TheDistanceToThePlane 


//  Think  of  a  tetrahedron  pd  at  the  peak,  P0,P1,P2  at  the  base 
//  if  the  Pd  projects  into  the  base  then  the  height  will  do. 

//  if  Pd  projects  outside  the  base  then  one  of  the  faces  will 
//  be  obtuse  to  the  base,  and  the  DistanceFromTheEdge  will  be 
//  required. 

if  (PdNormal*BaseNormal<0 . 0)  //  If  normals  are  not  in  the  same 

if  (d<Results . d)  //  direction  Or  if  the  distance 

{  f I  to  the  plane  is  more  then 

//  results  then  don't  bother. 


if  (BaseNormal*NormalToTriangle (PO, PI , Pd) <0 . 0) //  If  obtuse 

d=DistanceFromTheEdge ( Pd, PO , Pi ) ;  //  use  edge 

else  if  (BaseNormal *NormalToTr iangle ( PI ,  P2 ,  Pd)  <0 . 0) 
d^DistanceFromTheEdge (Pd,  PI ,  P2 )  ; 
else  if  (BaseNormal *NormalToTr iangle (P2, PO, Pd) <0 . 0) 
d=DistanceFromTheEdge (Pd, P2,  PO)  ; 

if  (d<Results . d)  //  if  distance  to  the  Triangle  is 

{  //  less  then  results  then  update  results 

Results. d  =  d; 

Results . Index=  Element . Index; 

} 

) 


Figure  1.  Pseudo  Code  for  Finding  the  Distance  to  a  Triangle 
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whether  or  not  the  triangle  Pj,  Po,  Pi  is  obtuse  at  Po,  or  P|,  if  it  is  then  the  distanee  from  the 
obtuse  point  to  Pj  is  the  distance  from  the  edge.  Otherwise  the  height  of  the  triangle  is  used  as 
the  distance  to  the  edge.  Figure  2  contains  pseudo  code  for  this  method. 

Interpolation 

Once  the  nearest  source  triangle  to  a  destination  point  has  been  identified  any  scalar  or 
vector  data  associated  with  the  points  of  the  source  triangle  can  be  interpolated  to  the  destination 
point.  This  application  uses  a  variant  of  baryccntric  coordinates  [3,4].  Baryccntric  coordinates 
arc  weighting  coefficients  that  when  multiplied  by  the  source  triangle  points  sum  to  the 
destination  point.  Specifically  baryccntric  coordinates  satisfy  the  following  equation, 

P^Co*Po+C,  *P,+C2*P2 
/=Co+C/+C2 

where  P  could  be  any  point  in  the  plane  of  the  triangle  and  Po,  Pi.  P2  arc  points  of  the  nearest 
source  triangle,  and  Co,  Ci,  and  C2  arc  the  baryccntric  coordinates.  By  definition  the  baryccntric 
coordinates  must  add  up  to  1 .  The  baryccntric  coordinates  arc  used  to  interpolate  any  scalar  by 
simply  multiplying  the  coordinates  by  their  respective  scalar  values  like  so, 

Sj=Co*S„+C,*S,+C2*S2 

Unfortunately,  for  this  application  Pd  is  not  likely  to  lie  on  the  plane  of  the  triangle,  so  it  will  be 
projected  to  the  nearest  point  in  the  plane  of  the  triangle.  The  baryccntric  coordinates  arc  then 
calculated  using  the  following  equations: 


Float64  BinaryTreeType : : 

DistanceFromTheEdge  (const  VectorTypeS  Pd, 

const  VectorTypeS  PO, 
const  VectorTypeS  PI) 

( 


VectorType 

VectorType 

VectorType 

Float64 


L  =P1-P0 
R0=Pd-P0 
Rl=Pd-Pl 
d; 


if  (R0*L<0.0) 

d=Abs (Pd-PO) ; 
else  if  (R1*L>0.0) 
d=Abs (Pd-Pl) ; 

else 

d=Abs (RO lUnit (L) ) ; 


// 

If 

triangle  obtuse 

at 

PO 

// 

Use 

distance  From 

The 

Pd 

to 

PO 

// 

If 

triangle  is  obtuse 

at 

PI 

// 

Use 

distance  From 

The 

Pd 

to 

PI 

// 

Use 

height  of  the 

tria 

ingle 

return  (d) ; 


); 


Figure  2.  Pseudo  Code  for  Finding  the  Distance  to  an  Edge 
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TotalArca=  AreaOfTriangIc  (Po,Pi,P2) 

Co  =  ArcaOfTrianglc(P,Pi,P2)/ TotalArca 
Cl  =  ArcaOfTriangle(P,P2,Po)/  TotalArca 
C2  =  ArcaOfTrianglc(P,Po,P|)/  TotalArca 


In  the  current  implementation  of  the  code  these  values  arc  always  positive  and  as  such  arc 
incorrect  if  P  is  outside  the  boundaries  of  the  source  triangle.  The  correct  coordinates  would 
need  to  have  at  least  one  negative  value  if  P  was  outside  the  boundaries  of  the  source  triangle. 
Using  the  baryccntric  coordinates  under  these  conditions  is  akin  to  extrapolation  rather  than 
interpolation.  In  order  to  avoid  this  type  of  extrapolation  and  to  ensure  that  we  don’t  use 
malformed  baryccntric  coordinates,  once  again  the  tetrahedron  (Pd,Po,Pi,P2)  is  considered,  and  if 
any  of  the  faces  arc  obtuse  then  the  method  is  switched  to  projecting  Pj  to  the  obtuse  edge  of  the 
base  triangle  and  calculating  the  coefficients  necessary  to  linearly  interpolate  Pj  from  the  edge 
points.  If,  however.  Pa  projects  outside  of  the  edge  points  then  the  nearest  point  in  the  triangle  is 
used.  Figure  3  contains  pseudo  code  for  calculating  the  coordinates. 


void  GeometryPartType : : 

UpdateBarycentricCoordinates (const  VectorType&  Pd, 

ElementType&  Element, 

SearchResultsType&  Results) 

{ 

Float64  C[3] ; 

Float64  TotalArea; 

VectorType  BaseNormal; 

VectorType  Pm; 

VectorType&  PO=Element . PO; 

VectorTypeS  Pl=Element . PI ; 

VectorType&  P2=Element . P2 ; 

BaseNormal=NormalToTriangle (PO,  PI,  P2)  ; 

Pm=Pd- ( ( Pd-PO ) *BaseNormal ) ^BaseNormal ;  //  Pd  projected  to  the  plain 

C [ 0 ) =AreaOfTriangle (PI , P2 , Pm) ; 

C[l] =AreaOfTriangle (P2, PO, Pm) ; 

C [2 ) =AreaOfTriangle (PO, PI, Pm) ; 

TotalArea  -  C [ 0 ) +C ( 1 ] +C [2 ] ; 

C[0)/=  TotalArea; 

C[l]/=  TotalArea; 

C[2]/=  TotalArea; 


Figure  3.  Pseudo  Code  for  Calculating  Barycentric  Coordinates 
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//  Think  of  a  tetrahedron  Pd  at  the  peak,  P0,P1,P2  at  the  base 
//  if  the  Pd  projects  into  the  base  then  the  height  will  do. 

//  if  Pd  projects  outside  the  base  then  one  of  the  faces  will  be 
//  obtuse  to  the  base,  and  the  DistanceFromTheEdge  will  be 
//  required. 

if  (BaseNormal*NormalToTriangle (PO , PI, Pd) <0 . 0 )  //  if  obtuse  face 

{ 

C[2]=0.0; 

UpdateBarycentricCoordinatesFromEdge ( Pd, PO , Pi, C [ 0 ] , C [ 1 ] ) ; 

) 

else  if  (BaseNormal*NormalToTriangle ( PI,  P2,  Pd) <0 . 0 ) 

{ 

C[0]=0.0; 

UpdateBarycentricCoordinatesFromEdge ( Po,  PI ,  P2 , C [ 1 ] , C [2 ] ) ; 

) 

else  if  (BaseNormal *NormalToTriangle ( P2 ,  PO ,  Po)  <0 . 0 ) 

( 

C[1]=0.0; 

UpdateBarycentr icCoordinatesFromEdge (Po, P2 , PO, C [ 2 ] , C [ 0 ] ) ; 

) 

Results . Barycentr icCoordinates [ 0 ] =C [ 0 ] ; 

Results . Barycentr icCoordinates [ 1 ] =C [ 1 ] ; 

Results . Barycentr icCoordinates [2 ] =C [ 2 ] ; 

]; 

void  GeometryPartType : : 

UpdateBarycentricCoordinatesFromEdge (const  VectorType&  Po, 

const  VectorType&  PO, 
const  VectorTypeS  PI, 

Float64&  CO, 

Float64&  Cl) 

i 

VectorType  L  =P1-P0; 

VectorType  R0=Po-P0; 

VectorType  Rl=Po-Pl; 

Float64  d; 

if  (R0*L<0.0) 

{ 

C0=1.0; 

C1=0.0; 

) 

else  if  (R1*L>0.0) 

{ 

C0=0. 0; 

Cl=1.0; 

} 

else 

( 

Cl= (R0*Unit (L) ) /Abs (L) ; 

C0=1.0-C1; 

) 

}; 


Figure  3.  (Continued) 
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Summary 

An  efficient  method  for  unstructured  mesh-to-mesh  interpolation  has  been  presented.  This 
method  uses  a  binary  space  partitioning  tree  to  sort  the  elements  of  the  source  mesh.  Using  this 
tree  data  structure  an  efficient  search  is  performed  for  the  source  mesh  element  nearest  each  of 
the  destination  mesh  points.  Once  found  the  nearest  source  mesh  element  is  used  to  compute 
baryccntric  coordinates  which  arc  then  used  as  weighting  coefficients  for  the  interpolation.  This 
method  has  been  successfully  implemented  and  used  for  coupling  of  a  CFD  and  CSD  code  for 
the  simulation  of  fluid  structure  interaction  problems. 
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